#ifndef HLLC_HPP
#define HLLC_HPP

#include "typeAliases.hpp"
#include "reconstruction_types.hpp"
#include "structMeshFunctions.hpp"
#include "macros.hpp"

/** @brief Harten, Lax, van Leer contact (HLLC) Riemann solver functions
 * 
 * Functions, whose parameters have the suffix L or R (for left/right) must be
 * provided with arguments pertaining to that specific side.
 * Functions with a WaveRegion template parameter produce different results
 * for left/right region, and must be provided with arguments matching the
 * specified WaveRegion.
 * Functions, whose parameters have a Side or OtherSide suffix produce the same
 * output when switching left/right arguments.
*/
namespace hms::hllc {

enum WaveRegion {
    left,
    right,
    star,
};


template<typename V, typename C>
auto vNormalStar(
	const V& vNormalL,
	const V& vNormalR,
	const C& cwL,
	const C& cwR
) -> decltype(auto);


template<typename V, typename C>
auto cwStar(
	const V& vNormalL,
	const V& vNormalR,
	const C& cwL,
	const C& cwR
) -> decltype(auto);


template<typename D>
auto pressureFlux( const D& d ) -> decltype(auto);


template<typename D, typename V>
auto massFlux(
	const D& d,
	const V& vNormal
) -> decltype(auto);


template<typename D, typename V>
auto momentumFlux_normal(
	const D& d,
	const V& vNormal
) -> decltype(auto);


template<
	typename S1, typename S2,
	typename F, typename Q>
auto hllFlux(
	const S1& sSide,
	const S2& sOther,
	const F& fSide,
	const F& fOther,
	const Q& qSide,
	const Q& qOther
) -> decltype(auto);


template<bool isOrtho, EdgeRange orien = EdgeRange::undefined>
HMS_INLINE
void setFlux(
	Ref<Array3s> fluxLeft,
	scalar sL, scalar sR,
	scalar dL, scalar dR,
	scalar vnL, scalar vnR,
	scalar vpL, scalar vpR,
	[[maybe_unused]] const Ref<const Vector2s>& normal
);


template<bool isOrtho, EdgeRange orien = EdgeRange::horizontal>
HMS_INLINE
void setFlux(
	Ref<Array3Xs> fluxLeft,
	const Ref<const Array1Xs>& sL,
	const Ref<const Array1Xs>& sR,
	const Ref<const Array1Xs>& dL,
	const Ref<const Array1Xs>& dR,
	const Ref<const Array1Xs>& vnL,
	const Ref<const Array1Xs>& vnR,
	const Ref<const Array1Xs>& vpL,
	const Ref<const Array1Xs>& vpR,
	[[maybe_unused]] const Ref<const Matrix2Xs>& normal
);


template<typename SolverType, EdgeRange orien = EdgeRange::undefined,
	typename F, typename Buf,
	typename D, typename VN, typename VP, typename N
	>
HMS_INLINE void hllcScheme(
	F& flux, Eigen::DenseBase<Buf>& buf,
	const D& dL, const D& dR,
	const VN& vnL, const VN& vnR,
	const VP& vpL, const VP& vpR,
	const N& normal
);

template<typename SolverType, EdgeRange orien = EdgeRange::undefined>
ArrayNNs<3,2> hllcScheme(
	const EdgeValues& ev,
	const ArrayNNs<2,2>& vRotLR,
	const Ref<const Vector2s>& normal
);


} // end namespace hms::hllc

#include "hllc.tpp"

#endif